#include <iostream>
#include <iomanip>
#include <cmath>
#include <complex>
#include <sys/time.h>
using namespace std;

const double pi = 3.14159265358979323846;
const int N = 128; // 多项式的最大支持位数
complex<double> temp[N]; // 用来充当临时数组

void generateSample(complex<double>* series) {
    for (int i = 0; i < N; i++) {
        int freq1 = 8;
        int freq2 = 16;
        int freq3 = 32;
        int freq4 = 64;
        series[i] = complex<double>(cos(2 * pi * i / freq1), sin(2 * pi * i / freq1));
        series[i] += complex<double>(cos(2 * pi * i / freq2), sin(2 * pi * i / freq2));
        series[i] += complex<double>(cos(2 * pi * i / freq3), sin(2 * pi * i / freq3));
        series[i] += complex<double>(cos(2 * pi * i / freq4), sin(2 * pi * i / freq4));
        //series[i] = (double)(rand()%1000000)/100000;// 随机生成一个0~10的小数
    }
}
void show(complex<double>* series) {
    for (int i = 0; i < N; i++) {
        cout << fixed << setprecision(5) << abs(series[i]) << ", ";
    }
    cout << endl;
}
void FFT(complex<double>* series, int len) {

    //递归出口:
    if (len == 1) {
        return;
    }
    //重新排列series:
    for (int i = 0; i < len / 2; i++) {
        temp[i] = series[i * 2];
        temp[i + len / 2] = series[i * 2 + 1];
    }
    for (int i = 0; i < len; i++) {
        series[i] = temp[i];
    }

    // 分治求A0和A1
    FFT(series, len / 2);
    FFT(series + len / 2, len / 2);

    // 通过A1和A2，计算A:

    double d_rad = 2 * pi / len; // 单位角幅度值
    for (int i = 0; i < len / 2; i++) {
        complex<double> x(cos(i * d_rad), -sin(i * d_rad));
        complex<double> temp1 = series[i];//A0
        complex<double> temp2 = x * series[i + len / 2];//A1
        series[i] = temp1 + temp2;
        series[i + len / 2] = temp1 - temp2;
    }
}
void DFT(complex<double>* series, int len) {
    for(int i = 0; i < len; i ++){
        complex<double> temp1(0, 0);
        for(int j = 0; j < len; j ++){
            temp1 += series[j] * exp(complex<double>(0, -2 * pi * i * j / len));
        }
        temp[i] = temp1;
    }
    for (int i = 0; i < len; i++) {
        series[i] = temp[i];
    }
}

int main() {

    //产生随机时序信号序列:
    complex<double> a[N];
    complex<double> b[N];
    generateSample(a);
    for (int i = 0; i < N; i++) {
        b[i] = a[i];
    }
    cout << "时序信号：" << endl;
    show(a);

    //计时:
    struct timeval start1;
    struct timeval end1;
    struct timeval start2;
    struct timeval end2;

    //FFT:
    gettimeofday(&start1, NULL);
    FFT(a, N);
    gettimeofday(&end1, NULL);
    //DFT:
    gettimeofday(&start2, NULL);
    DFT(b, N);
    gettimeofday(&end2, NULL);

    //用时统计:
    double time1 = (end1.tv_sec - start1.tv_sec) + double((end1.tv_usec - start1.tv_usec)) / 1000000;//单位s;
    double time2 = (end2.tv_sec - start2.tv_sec) + double((end2.tv_usec - start2.tv_usec)) / 1000000;//单位s;
    cout << "FFT用时：" << time1 << "s" << endl;
    cout << "DFT用时：" << time2 << "s" << endl;

    //输出结果:
    cout << "FFT结果：" << endl;
    show(a);
    cout << "DFT结果：" << endl;
    show(b);
    return 0;
}
